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I. INTRODUCTION 

The charge dynamics during the bombardment of di¬ 
electric samples by the primary electrons (PE) in a 
Scanning Electron Microscope (SEM) can be studied 
with diff erent met hods. One is the Monte Carlo (MC) 
techniqu^^^i^IloISnHlIl^ whose main advantages are the rig¬ 
orous semi-classical account of the microscopic physics 
and the ability to consider non-equilibrium dynamics (by 
this we mean the dynamics of particles with their energies 
not yet distributed as in thermal equilibrium). However, 
the MC approach is known to suffer from the increase of 
computational complexity in the case of long-range po¬ 
tentials, such as those of the electrostatic fielcP. More¬ 
over, in an inhomogeneous sample the required potentials 
can only be obtained by numerically solving a large elec¬ 
trostatic problem with a very high and non-uniform spa¬ 
tial resolution at each step of the MC algorithm. Also, for 
a reliable estimate of the particle flux through a part of 
the sample boundary, one needs to consider a sufficiently 
large statistical ensemble, which is computationally ex¬ 
pensive. 

For these and other reasons an alternative and in many 
ways a much simpler self-consistent approac h originating 
in semiconductor physics has been proposecP^^tlllMHl 
This so-called Rostoc Program takes the current density 
point of view, considering currents rather than charge 
densities to be the fundamental quantities. Some of 
the advantages of the current-based approach are: the 
possibility to model the sample-vacuum interface via a 
reflection-transmission coefficient formalism and to in¬ 
clude the tertiary electrons returning to the sample into 
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the model. On the other hand, it is more difficult to de¬ 
scribe proper ohmic contacts in this way and it is hard 
to extend this approach to two and three spatial dimen¬ 
sions. 

The traditional self-consistent approach of the semi¬ 
conductor physic^ considers the charge densities obey¬ 
ing the drift-diffusion-reaction (DDR) system of equa¬ 
tions to be the fundamental quantities. This approach is 
particularly suited for modeling the equilibrium charge 
transport in inhomogeneous semiconductor devices (e.g. 
junctions). Some parts of the D DR model ha ve already 
been applied to the SEM problenPISMllIlIIll] However, 
these previous studies have omitted the trapping rate 
equations thereby missing an important feature of the 
charging dynamics. Also, the model employed relies on 
an MC treatment of the primary electrons (PE), their 
initial scattering, and the emission of the secondary elec¬ 
trons (SE) through the sample-vacuum interface. Hence, 
the question remains whether a fully self-consistent con¬ 
tinuum DDR model without any MC parts can ade¬ 
quately describe the charging of dielectric samples by a 
focused electron beam. 

The main challenges one faces in developing a fully 
self-consistent DDR model for the SEM problem are: the 
non-equilibrium charge injection mechanism followed by 
the generation of secondary particles via ionization, the 
fact that secondary electrons may leave via the vacuum- 
sample interface, the back-coupling effect of the accumu¬ 
lated charges on the primary beam, and the multi-scale 
nature of the problem (spatial as well as temporal). Here 
we show that all these problems can be successfully solved 
and that the traditional DDR approach represents a vi¬ 
able alternative to MC simulations. 

There are obvious limitations to the classical equilib¬ 
rium continuum picture of the particle dynamics inside 
dielectrics. For example, one does not expect the DDR 
approach to be applicable on the level of a single PE or 
in the first moments following its impact. Yet, the MC 
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simulations indicate that the cloud of secondary parti¬ 
cles created by ionization contains sufficiently many par¬ 
ticles so that their subsequent evolution can indeed be 
described on the level of densities. Moreover, the MC 
simulations and the many controlled experiments pro¬ 
vide sufficient information to construct a semi-empirical 
source function that mimics the impact and its immedi¬ 
ate aftermath for each PhP. 

As far as the exit of SB’s through the sample-vacuum 
interface, the concept of the surface recombination ve¬ 
locity (SRV) appears to be sufhciently robust to describe 
this proces^. This SRV, roughly speaking, determines 
the rate at which particles are allowed to leave and de¬ 
pends on the materials adjacent to the interface. Due to 
the virtual absence of data for dielectric-vacuum inter¬ 
faces, however, the SRV remains a tuning parameter in 
our method. 

In our approach, a set of equations is employed for 
the recombinati on and trapp ing rates, whereas, previ¬ 
ous DDR studie^^^^mHIllIM gjjjiply use fixed values for 
these rates. We show that the trapping of particles intro¬ 
duces a large-scale (slow) dynamics into the picture and 
determines not only the main features of the charge den¬ 
sity distribution inside the sample, but also the abrupt 
changes in the surface charge density prior to the estab¬ 
lishment of the steady state. 

The self-consistent DDR approach presented here 
brings its own set of unique challenges and questions with 
it. For example, one has to take care that the numerical 
solver is sufficiently robust and stable and does not pro¬ 
duce non-physical (e.g. negative) values for the particle 
concentrations. Also, a purely theoretical question arises 
about the existence of a steady-state and/or periodic so¬ 
lutions to the DDR equations. 

The remainder of this paper is organized as follows. In 
the next section we describe the set of equations pertain¬ 
ing to the drift-diffusion model. Then, a separate section 
is devoted to the model of the charge injection process. 
A section and an Appendix describe the details of the nu¬ 
merical solution via the Finite-Element Method. Finally, 
a series of numerical experiments is presented followed 
by the conclusions. 


equation: 


■ {e\/V) =-{p + pt-n-nt), (1) 

£0 

where q is the elementary charge, n(x, t) is the density of 
free electrons, nt{x,t) is the density of trapped electrons, 
p(x, t) is the density of free holes, pt(x, t) is the density of 
trapped holes, the constant Eq is the dielectric constant 
of vacuum, and the function e(x) is the (static) relative 
permittivity of the sample. 

The continuity and trap rate equations can be stated 
as 


dn ^ ^ rr ^ 

yw + ^- Jn = U + Sn 
at 


dnt 
dt ’ 


( 2 ) 


dnt 
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( 3 ) 
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U + Sp 


dpt 
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( 4 ) 


dp 

= fJpVthiNp - pt){p - n^) - -fppt, (5) 

with the constitutive relations for the current densities 
given by 


Jn — Dji^n -f /iT^uW, 


( 6 ) 


Jp — DpS/p |XppW, 


( 7 ) 


where Pn ^^nd Pp are the electron and hole mobilities, 
Dn and Dp are the diffusion coefficients, cr„ and Up are 
the electron and hole trapping cross sections, 7 „ and 7 p 
are the detrapping time constants, and Np are the 
densities of trapping sites, Vth is the thermal velocity, 
Ui is the intrinsic carrier density, and Sp are source 
functions which will be defined in section |III A and U 
is the charge recombination rate given by formula (|^ in 
section III Cl 


II. DRIFT-DIFFUSION-REACTION MODEL 
A. Basic equations 

The DDR model consists of a set of three coupled non¬ 
linear PDFs and two nonlinear ODEs. Namely, the po¬ 
tential equation, two continuity equations (one for the 
electron and one for the hole current densities), and two 
trapping rate equations for trapped electrons and hole^^. 
Thus, we monitor the simultaneous space-time evolution 
of four species of particles and one potential function. 

The electrostatic potential V (x, t) satisfies the Poisson 


B. Trapping and detrapping 

The process that causes low-energy charges in di¬ 
electrics to be transferred to a localized state is called 
trapping. Trapping occurs at a trapping site. The 
charges that have been trapped at a certain site at 
one time, due to several reasons, for instance, the field- 
induced detrapping, can get detrapped and become free 
at a later time. The process can continue which means, 
this free charge can get trapped again somewhere els^i^. 
A detailed analysis of the electron and hole trapping in 
dielectrics can be found iiP. In the present model, this 
process is described by the two ordinary differential equa¬ 
tions ^ and ([^. The coefficients UnVth (o'pVth) and 7 „ 
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(7p) specify the rate of electron (hole) trapping and de¬ 
trapping, respectively. It is easy to foresee that initially 
the terms drit/dt and dpt/dt in equations (§ and 0 
will act as time-dependent sink terms. However, as soon 
as the density of trapped charges reaches the density of 
trapping sites {Nn and Np) or the density of free particles 
drops below n^, these terms will act as time-dependent 
sources. 


C. Charge recombination 

There are two basic recombination mechanism in 
semiconductor physics described by the Auger and the 
Shockley-Read-Hall (SRH) model^^. It is known that 
the Auger model is more appropriate at higher carrier 
concentrations caused, e.g., by heavy doping or high level 
injection under concentrated sunlight. Therefore in the 
present case, where the concentrations are not that high, 
we opt for the SRH model. 

The function U(n,p) in ([^ and Q is the generation- 
recombination rate, in other words, the rate at which 
electron-hole pairs are generated minus the rate at which 
they are recombined. Since electrons and holes are gen¬ 
erated and recombined in pairs, we have the same rate 
function for the two species. In the SRH modeP^ this 
function is given by 


U{n,p) 


nf — np 

T„(n -I- Ui) + Tp{p + m) ’ 


( 8 ) 


where t„ and Tp are the life time parameters for the elec¬ 
trons and holes, respectively. 


D. Boundary and initial conditions 

The SEM chamber consists of two main parts - the 
vacuum and the sample. Considering a cross-section, 
we assume a rectangular outer boundary Fig. which 
can be further adjusted to take the actual geometry into 
account. The domain is further divided into two equal 
parts, where one represents the sample and the other the 
vacuum chamber. The Poisson equation 0 is considered 
on the whole domain (Hi and H2), whereas, equations 
0-0 are solved on the lower domain (f22) only. 

Depending on the material in contact with the sample 
two types of boundary conditions are common: Dirich- 
let conditions at ohmic contacts, and Robin conditions 
at Schottky and similar semi-insulating contactd^. The 
boundary of the sample domain D2 consists of the Dirich- 
let part 9D2 (where the sample is in contact with the 
walls of the SEM chamber or another highly conducting 
material), and the Robin part S (sample-vacuum inter¬ 
face). At ohmic contacts (sides and the bottom of D2) 
the space charge vanishes, i.e., 

p— n = 0 on (9D2 X [0,tend]- (9) 


Furthermore, the system is in thermal equilibrium there, 
which is expressed by the relation 

np = nj on (9D2 x [0,tend]- (10) 

From the above relations, we have 

n(x, t) = m, _p(x, t) = m on 5D2 X [0,tend], (H) 

We also assume homogeneous Dirichlet conditions for the 
potential on the wall of the SEM chamber, i.e. 

F(x,t) = 0 on (5Di U 9D2) X [0,tend], (12) 

which could be easily adjusted to account for any finite 
value of the electric potential. 

The dielectric-vacuum surface recombination model 
can been obtained as a simplification of the SRH modeP^ 
and leads to a somewhat unusual Robin-type boundary 
condition at the sample-vacuum interface. Namely, it is 
a semi-insulating contact for the electrons and an insu¬ 
lating contact for the holes (since holes cannot exist in 
vacuum): 


*In ‘ ^ ^i) OCL S X [0,tend]; (^8) 

Jp . 1/= 0 on Sx[0,tend], (14) 

where is the surface recombination velocity (SRV) for 
electrons and v denotes the unit outward normal vector 
on the boundary E. This parameter has an important 
role in the model and will be discussed later. In fact, the 
insulating condition ( |l4| ) can be also considered a semi- 
insulating contact with its surface recombination velocity 
set to zero {vp = 0). 

The intrinsic carrier density has been considered as the 
initial condition 

n(x,0)=ni, p(x, 0) = nj, in D2. (15) 


III. BEAM MODEL 

A. Impact of an individual primary electron 

When an electron beam illuminates the sample some 
of the primary electrons will reflect as backscattered elec¬ 
trons (for silicon oxide on average 20%), while the rest 
penetrates the sample and produces a large number of 
secondary electrons/holes. 

It is important to realize that the DDR model Q-(^ 
describes the equilibrium transport of charged particles 
with the distribution of kinetic energies either depending 
only on the (effective) temperature of the sample or sim¬ 
ply being stable throughout the simulation time. How¬ 
ever, the impact of the primary electron and its immedi¬ 
ate aftermath are not equilibrium processes. Yet, since 
typical dielectric samples are made of dense materials, 
the numerous collisions will lead to the thermalization 
of all generated secondary particles shortly after the PE 
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function: 
L(t) = 




— kt ’ 


dL 

dt 


= kL{l-L). (17) 


where k is the Malthusian parameter and w is an ini¬ 
tial condition related to the so-called carrying capacity 
ranging from 0 to 1. The values for k and w used in our 
calculations are reported in Table |T] We choose L to be 
the logistic function since pair creation is an avalanche- 
type process and as such is mathematically similar to the 
population growth. 

In (16) tg denotes the thermalization time, which is 
taken here to be approximately the time of the ballistic 
flight of the primary electron. Special relativity provides 
a simple relation between the velocity of a primary elec¬ 
tron and its energy: 


V = cjl- 


1 


(1 + 


(18) 


FIG. 1. General schematics of the problem. 


impact, so that their subsequent transport can indeed be 
modelled by §" 0 - The precise rate of this thermaliza¬ 
tion is not known, but could be obtained with dedicated 
Monte-Carlo simulations^, which are beyond the scope 
of the present paper. Here, we simply assume that all 
the secondary particles are in equilibrium by the time 
the primary electron looses most of its kinetic energy. 

To account for the initial non-equilibrium transport we 
introduce an effective source model based on the avail¬ 
able experimental data about the secondary charge dis¬ 
tribution inside dielectrics. Moreover, since most of the 
SE emission happens during this initial non-equilibrium 
stage (in practice as well as in our simulations), the sur¬ 
face recombination velocity of the sample-vacuum in¬ 
terface (discussed below) allows further fine-tuning of 
the model using the material-dependent experimental SE 
yield data. 

Mathematically, the injection of electrons is described 
by the terms S'„(x, t) and S'p(x, t) in the right hand sides 
of the continuity equations ([^ and Q. In our model 
these source functions consist of two factors. The first 
factor represents the density of charge at the end of the 
ionization stage following the impact by a primary elec¬ 
tron. The second purely temporal factor approximates 
the dynamics of the ionization stage, i.e., the build-up 
of the secondary charge during the first picosecond after 
collision. Thus, the source function has the form 


Sn,p (^7 t) 


^Ti ,p ) dL 

L{tg) — L{0) dt ’ 

0 , 


if 0 <t <tg; 
otherwise; 


(16) 


where ft.(x. Egg) is the charge distribution function de¬ 
pending on the effective energy of the primary electron, 
as will be explained shortly, and L is the following logistic 


where c is the speed of light in vacuum. The time of 
flight tg can be estimated by dividing the penetration 
depth (will be explained below) by this velocity (or a 
twice lower ‘average’ velocity). In either case it appears 
that for the relevant range of primary energies tg is in 
the order of 10“^^ seconds, i.e., extremely short with 
respect to the average time between electron impacts in 
a typical SEM beam. If this estimate is correct, then the 
DDR model is indeed applicable to the charge dynamics 
not only on large time scales, but also on the scale of 
individual impacts. 

Let R{Eq) denote the maximum penetration depth by 
the primary electrons with initial energy Eq. There exist 
several empirical formulas for R(Eq). For example, the 
experimental results by Pott^^ indicate that R is given 
by: 

^1.5 

R = 0.1 ° [itm], (19) 

P 


where p is the mass density of the sample material in 
g/cm^ and Eq is in keV. On the other hand, theoretical 
considerations by Kanaya and OkayamgP^ lead to 

4 ^ 1.67 

R = 2.1QxlQ-^-^[pvnl (20) 

where Z is the atomic number, A is the atomic mass and 
Eq is in keV. The following composite form ula propose d 
by Fitting has been used by several author 

/ 900p-O-8£;i-3 [A] for Eq < 10 keV, 

I [A] for Eq > 10 keV, ^ ’ 


where p is in g/cm^. Here we employ the most recent 
estimate by Fittin^^il: 


^1.45 

i?(p,Eo)=93.4^ [nm], (22) 















5 


where p is in g/cm^ and Eq is in keV. 

According to several studie#*S*24| ^;]^e actual distri¬ 
bution of the secondary electrons and holes is well- 
approximated by a three-dimensional Gaussian function 
with its focus Xq located 0.3i? below the vacuum-sample 
interface: 


5(x, Eo) 


aA 


exp {—B\x 


xqI 


(23) 


additional particle upon the integration (28). Thus, the 
factors of eq. (16) are 


A' 0 ff ) 


^11.58 


EeS 


13.158 


X 


1 

ttR^ 


exp 




A-eff) — A-eff), 


(29) 


where Ei is the mean creation energy for one SE, a is the 
yield factor close to one, and 


7.5 


i? 2 ’ 



(24) 


where Egg depends on the surface electric potential at 
the time of impact. Of course, the source functions pro¬ 
posed here are only approximations. Nevertheless, they 
are based on the best experimental evidence and first 
principles calculations available to date. 


The constant C{Eq) is proportional to the fraction 77 of 
backscattered PE. For silicon, silicon dioxide, and alu¬ 
minium oxide, with 77 « 0.2, C can be obtained from: 


B. Bombardment and temporal smoothing 


C = 1.544^, (25) 

K 

where C is in eVA ^ and Eq is in keV. 

To account for the action of the surface potential I 4 on 
the primary electron, we introduce the effective energy 
Egg = Eq + Vs{ti), where ti is the time of impact, which 
should be applied in the distribution function instead of 
Eq, thus arriving at: 


g(pc,Egg) = 11.58 


Egg 

nR^Ei 


7.5, , 

I “ ^o| 


i? 2 ' 


(26) 


where R{p,Egg) is given by (22) and the pair creation 
energy E, depends on the material of the sample viaP^ 


Ei Ki 3 Eg + 1 eV, 


(27) 


with Eg denoting the energy gap of the material in eV. 

The total numbers Nse,sh of secondary electrons and 
holes corresponding to the distribution (26) can now be 
estimated as 


Nse = Nsh 


B3_z>0 


E ff 

5 (x,Eeff)dE« 0.877^, ( 28 ) 


showing that approximately 88 % of the effective energy 
is spent on the creation of charge pairs, which gener¬ 
ally agrees with MC simulations. According to (28) the 
number of secondary electrons generated by one primary 
electron is somewhere between hundreds and thousands. 
Hence, we may expect the drift-diffusion-reaction ap¬ 
proach to be a reasonable approximation at this scale. 

Thus we take the hp for holes in (16) to be equal to g as 
introduced in (26). Whereas for the electrons we recall 
that the primary electron is still present in the sample 
at t = tg. Hence, we adjust the coefficient in front of 
the exponent in the function g so that it features one 


Depending on the beam current primary electrons may 
arrive at an average rate as high as tens of millions 
per second. Previous applications of the drift-diffusion- 
reaction approach typically describe the SEM beam as 
a constant flux of electrons. The goal of the present 
paper is to avoid the latter approximation and directly 
consider, say, m, primary electrons arriving at times U, 
i = 0,1,... , 777 . Thus, one obtains a pulsed source where 
the next PE arrives in a medium with some residual 
charge left from the impact of the previous PE. 

Although, we gain some valuable insights about the 
subsurface charge dynamics and the effect of beam cur¬ 
rent, it is obviously too time consuming to consider bom¬ 
bardments of a sample by a large number of electrons in 
this way. Hence, a different approach is needed to study 
saturation effects at larger time scales. Also, the SE yield 
calculations on the level of single PE’s, although possible, 
are hard to justify and interpret. 

The main technical challenge preventing direct large- 
scale simulations with our method is the pulsed nature 
of the source terms requiring many time steps to be per¬ 
formed by the solver between electron impacts. A way to 
reduce the computational burden is to derive a smoother 
function describing the behavior of source terms at larger 
time scales. In the limit such a smoother source function 
should approach the constant beam currents of the other 
DDR models. 

To achieve this we employ a temporal average of our 
source function, which also mimics the way the SEM 
response is measured (time-averaged yield, rather than 
the yield due to individual PE’s). The average value of 
S'„(x, t) over a period of time T between the impacts can 
be expressed as 


5'n(x) = [ Sn{x,t')dt', (30) 

^ Jo 

and is a time-independent function. In what follows we 
call this a time-uniform or simply a uniform source. 
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Unfortunately, smoothing of the source has its price. 
Due to the presence of nonlinear terms in Q-0, so¬ 
lutions obtained with a time-averaged source term will 
not be the exact time-averaged values of the unknowns, 
but only the approximations thereof. Hence, to apply the 
DDR approach at both time scales successfully one needs 
to define constitutive relations and material parameters, 
such as the surface recombination velocity, for each scale 
separately. This is the so-called homogenization problem, 
typical for spatial multiscale analysis in physics (e.g. ef¬ 
fective medium problem in electrodynamics). 

Further, although a uniform source switched on at 
t = 0 may be expected to eventually produce a steady- 
state distribution of charge, it is an open theoretical ques¬ 
tion whether the actual pulsed source leads to the corre¬ 
sponding periodic charge variations. 


IV. NUMERICAL METHOD 
A. Numerical scaling 


Z 



To avoid numerical difficulties and maintain the accu¬ 
racy of the solution, a simple scaling of variables has been 
performed. To this end we introduce a set of characteris¬ 
tic dimensionless quantities. We denote the characteris¬ 
tic length scale by /*, the characteristic time scale by t* 
and the characteristic density scale by p*. The numerical 
values of these dimensionless parameters relevant to the 
scale of the present problem are 

r = 10"®, t* = 10"^^ p* = 10^®. (31) 


There is a relation between these values {t* = (U)^ and 
p* = (U)"®) that doesn’t change the form of the equa¬ 
tions, so that one only needs to introduce the rescaled 
versions for some of the constitutive parameters: 




‘^iP 


rii = —, e = 


t* ’ p*’ - p*{i*Y' 

^ P ^n,pi ^n,p — ^ '~fn,p') 


iV - 

^ n..r> — 1 


S=—S= ^(?(i), 
p* dt ’ 


Also the boundary and initial conditions should be 
rescaled, since, e.g. the rescaled version of the surface 
recombination velocity is given by: 


Vn 


l*V„ 


(33) 


B. FEM solver 

We employ the finite element method (FEM) for the 
numerical solution of the coupled system 0-0 and im¬ 
plement it as a solver within the COMSOL Mu tiphysics 
package. Although, there are many predefined modules 
and solvers in COMSOL, none of them can be directly 


FIG. 2. Cylindrical geometry. 


applied with the present problem. The closest match is 
the semiconductor module. However, it is neither suited 
for studying the two different domains defined above, i.e., 
Oi U O 2 for equation 0 and fl 2 for the rest, nor does 
it allow to incorporate the additional equations (|^ and 
([^. Therefore, we have opted for building a new model 
using the general PDE and the ODE/DAE interfaces of 
COMSOL. 

Since the charge densities may be extremely concen¬ 
trated around the impact zone and form very thin layers 
near the vacuum-sample interface, a careful discretiza¬ 
tion strategy is required. To achieve sufficient accuracy 
one is advised to use the adaptive mesh with refinement 
in the impact zone and at the interface as well as the 
second-order Lagrange shape functions. A fully cou¬ 
pled approach with Newton-Raphson solver and adaptive 
time-stepping algorithm has shown the best performance 
with the present problem. 

The computational complexity of the problem pro¬ 
hibits a full three-dimensional (3D) simulation of realistic 
domains with sufficient spatial resolution on a standard 
PC. Nevertheless, one can obtain 3D results for certain 
configurations at a typical two-dimensional (2D) cost by 
exploiting their symmetry. Consider, for example, the 
cylindrical geometry presented in Figure In the cylin¬ 
drical coordinate system (r, 9, z) the PE beam imping¬ 
ing along the z-axis corresponds to the source term and 
boundary conditions independent of the angular coordi¬ 
nate 6. The solution will also be independent of 9 and 
the original 3D model is reduced to a 2D model in the 
(r, z)-coordinates. To arrive at the corresponding FEM 
solver the PDE’s 0,0 and 0 must be written in the 
so-called weak form, which is derived in the Appendix. 
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TABLE I. Parameters of dielectric materials. 


Parameter 

SiOa 

AI 2 O 3 

Unit 

e 

3.9 (ref.[^ 

10 (ref.l22) 



20 (ref.l^ 

4 (ref.li2) 

cmW-"s-^ 

Up 

0.01 (ref.l^ 

0.002* 

cmW-"s-i 

^71 

10"^® fref.l^ 

10-^® 

cm^ 

(T p 

10"^® (ref.l3Sl) 

10-18 

cm^ 

Vth 

10'^ (ref.l^ 

10" 

cm s”" 

'^n,p 

2 X 10“® fref.ESl) 

2 X 10“® 

s 

p 

2.65 (ref.l^ 

3.98 (ref.l22) 

g cm”® 

Es 

9 (ref. 122) 

9 

eV 

Nn,p 

1.6 X 10^® fref. 132) 

1.6 X 10^® 

cm”® 

yn.p 

10"* (ref.li2) 

10"^ 

S-" 

k 

25 

25 

S-" 

w 

10“® 

10“® 



*This value could not be found in literature and has been 
chosen by analogy with the relation between the electron 
and hole mobilities in Si 02 . 


V. NUMERICAL EXPERIMENTS 


fact that it is one of the few directly measurable quanti¬ 
ties in SEIVPS. Monte Ca rlo simulations of the SE yield 
are also availabl^^^^SIiniSI]^ Therefore, in the tuning of the 
SRV parameter one would mostly be relying on the SE 
yield data as a function of the PE energy. In the present 
continuous approximation the SE yield is computed as 
the flux density through the boundary E integrated over 
this boundary and over time from t = 0 to t = tend, and 
divided by the number of PE’s that arrived at the sample 
during that time interval. 

The following steps describe a simple optimization pro¬ 
cedure for tuning the value of SRV: 

• Let Vixp be the SE yield measured fro PE’s with 
energy Eq. 

• Let be the initial guess for the SRV, and let 
Y(vn'^) be the SE yield computed by the DDR 
solver with the SRV set to . 

• For Vn^ sufficiently close to the true (optimal) value 
we can assume a linear relation: 


A. Parameters and parameter-tuning 


It is clear that the values of the many constitutive 
parameters in the equations 0-0 may have consider¬ 
able influence on the results of simulations. In the ideal 
situation these parameters are either measured in dedi¬ 
cated experiments or computed from the first principles 
of quantum physics. While this status-quo has long been 
established with the usual semiconductor materials, the 
data for insulators are virtually absent. Often, very dif¬ 
ferent values for the same parameter are reported in the 
literature, which could be the result of varying sample 
properties, experimental conditions, or even trivial hu¬ 
man errors (see e.g. detrapping rates iiP^ anJMl) jn 
particular, there is a lot of uncertainty about the pa¬ 
rameters of recombination processes in the bulk and at 
interfaces. With this in mind we have made a selection 
of typical values for two distinct dielectrics shown in Ta¬ 
ble H 


In the boundary condition (13) the surface recombi¬ 
nation velocity Vn plays an important role analogous to 
the reflection coefficient of the current-based approach. 
In the literature the SRV is mostly discussed and mea¬ 
sured for the metal-oxide-semiconductor interfaces. With 
respect to our case, which is a vacuum-semiconductor 
interface, the only reference that could be found refers 
to Germanium under ion bombardmenlP^, reporting the 
value of (5 to 7) x 10^ cm/s. Considering the SRV to be 
material-dependent we believe that its value should be 
determined based on the consistency between the results 
of the present model and dedicated experiments. 

The SE yield is defined as the number of secondary 
electrons emitted through the sample-vacuum interface 
and picked up by the detector per one incident elec¬ 
tron. The importance of the SE yield stems from the 


Vexp - Y =a(v„- (34) 

Since, obviously, Y (0) = 0, the coefficient a can be 
obtained as a = ^ , so that Vn^ = Vn^ wgv ■ 

Vn ’ J 

• Compute Y(vn^) with the DDR solver. 

• If V(ui^^) is sufficiently close to Vexp, then stop and 

set Vn = Vn'^. Otherwise, continue with Vn'^ as the 
new initial guess. 

In principle, this process should be repeated with the SE 
yield data for a whole range of PE energies Eq. Unless 
changes in Eq significantly alter the temperature of the 
sample, the SRV of a given material is supposed to be 
independent of Eq. Hence, if the corresponding tuned 
values of for some material are all close to each other, 
then we have an additional confirmation that the DDR 
method is working properly. 

We conclude this section with a word of caution con¬ 
cerning the use of SE yield in determining the SRV. First 
of all, typical SE yield data correspond to some kind of 
stationary regime. It is known, however, that the SE 
yield keeps changing after the start of bombardment for 
quite a long time. Hence, we can compare the results 
of simulations with experimental data and tune the 
parameter only upon bombardment of the sample with 
a sufficiently large number of PE’s. In our approach 
such long-time simulations can only be performed with 
the time-uniform source. This means, however, that the 
tuned value of will be effective in nature. 

Secondly, numerical experiments demonstrate that the 
electron flux through the sample-vacuum interface is not 
only time-dependent, but also depends on the extent of 
the sample and the proximity of ohmic contacts or other 
conducting materials to the beam’s entry point. Hence, 
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FIG. 3. Evolution of the total charge {p + pt n nt)q FIG. 4. Evolution of the total charge {p + pt — n — nt)q 
(G cm-®) in Si 02 after the impact of a single primary electron (c cm"®) in AI 2 O 3 after the impact of a single primary elec- 


with the energy Eq = 1 keV. Top row: vertical cross-section, 
side length is 100 nm; bottom row: top-view of the sample 
vacuum interface, diameter is 200 nm. 


tron with the energy Eq = 1 keV. Top row: vertical cross- 
section, side length is 100 nm; bottom row: top-view of the 
sample-vacuum interface, diameter is 200 nm. 


one may expect different SE yield values with different 
samples of the same material. 


B. Impact of a single primary electron 

In this section we investigate the events following the 
injection of a single primary electron into a neutral di¬ 
electric sample. The goal of these numerical experiments 
is to estimate the space-time scales of the dynamics sepa¬ 
rately for all four particle species, i.e., n, p, nt, and pt, as 
well as the total charge density q{p + pt — n — rit) and the 
potential V. In particular, these simulations will help 
us to demonstrate that despite the poor mobility and 
diffusivity of dielectrics, the drift and diffusion of free 
charges is generally much faster than the characteristic 
time scales of the charging process. 

We focus on two common materials featured in many 
of the previous studies, namely, on the oxides Si 02 and 
AI 2 O 3 . As one can see from the data of Table |l] the 
difference between these materials is in the values of the 
relative permittivity e (Si 02 has a smaller e), the electron 
and hole mobilities (Si 02 is relatively more conduc¬ 
tive), and the mass density p (Si 02 has a smaller density). 
The effective values of the SRV, = 100 cm/s for Si 02 
and Vn = 200 cm/s for AI 2 O 3 , were obtained with the op¬ 
timization procedur e expla ined above and typical experi¬ 
mental SE yield datcP^^I^ using the time-uniform source 
model. Below we focus on a fixed PE energy Eq = 1 keV. 

Figures and 1^ show the snapshots of the time evolu¬ 
tion of the four charge species in the two materials. As 
can be seen from the upper-left images of these figures, 
the smaller mass density of Si 02 means, see eq. ( 22 ), that 


tial charge distribution are different as well, see eq’s. (26) 


with the same Eq the maximum PE penetration depth R 
and the center of the initial charge distribution are deeper 
for Si 02 than for AI 2 O 3 . The overall shapes of the inl¬ 


and (29), with the one of Si 02 being broader. Hence the 


DDR dynamics starts with different initial states in these 
materials. 

The generation of charge pairs by ionization takes place 
in the period of 1 picosecond after injection (tg=l ps). 
Note that in our model the processes of recombination 
and trapping begin already at t = 0. At tg the den¬ 
sity of free electrons is already beginning to decrease. In 
fact, the density of free electrons reaches its maximum 
of 2.07 X 10^® cm”® and 6.35 x 10®® cm”® for Si 02 and 
AI 2 O 3 , respectively, at around t = 0.6 ps. The density of 
free holes reaches its maximum roughly at t = 0.7 ps and 
remains constant till the end of the generation period tg. 
The maximum density of free holes in this stage for Si 02 
is 2.31 X 10®® cm”® and for AI 2 O 3 is 7.01 x 10®® cm”®. 

As expected, the density of trapped electrons and holes 
initially increases with time reaching, respectively, (at 
tg) the values of 1.79 x 10®®’ cm”® and 2.02 x 10®^ cm”® 


in Si 02 , and 5.31 x 10®® cm”® and 6.09 x 10®® cm”® m 
AI 2 O 3 . The lower density for the trapped holes compared 
with trapped electrons is due to the smaller trapping 
cross sections. In Si 02 the density of trapped electrons 
reaches its maximum of 2.01 x 10®® cm”® at f = 50 ps. 
Then, for more that 50 ns, which is a relatively long 
time, no change is seen in the distribution of trapped elec¬ 
trons. After that the maximum density of trapped elec¬ 
trons starts to decrease dropping to 1.69 x 10®® cm”® at 
t = 1 /is. For the trapped holes, reaching the maximum 
density takes much longer time compared to the trapped 
electrons. The density of trapped holes keeps increasing 
at f = 1 /xs and reaches its maximum of 1.71 x 10®® cm”® 
at t = 50 ns. The maximum density for trapped holes at 
time t = 1 ns is 2.74 x 10®® cm”® and at t = 1 /xs the 
density of the trapped holes is 1.7 x 10®® cm”®. A similar 
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dynamics of trapped particles is observed in AI 2 O 3 . 

During the first microsecond a fast decrease in the den¬ 
sity of free electrons and a slower decrease in the density 
of free holes is observed owing to the relatively strong 
trapping of electrons and a weaker trapping of holes. 
In Si 02 the major drop in the density of free electrons 
happens during the first 50 ps. The maximum density 
of free electrons is 5.21 x 10^'^ cm“^ at t = 50 ps and 
reaches almost the intrinsic carrier density of the mate¬ 
rial at t = 1 /rs. At t = 50 ps the density of free holes is 
higher than that of the free electrons ( 2.11 x 10 ^® cm“®). 

The interplay of the four charge species leads to the 
total charge density resembling an expanding spherical 
wave with initially a positive charge region in the mid¬ 
dle surrounded by a shell of negative charge followed by 
a negative middle region with a positive shell. Due to 
the emission of electrons the initial predominantly nega¬ 
tive charge at the sample-vacuum interface is gradually 
replaced by the positive charge. In Si 02 , the positive 
charge reaches its maximum of 0.13 C/cm® at the surface 
at t = 50 ps and the negative charge has the maximum of 
(in the sens of absolute value) —0.05 C/cm® at t = 2 ns 
and is situated close to the surface . In AI 2 O 3 , the pos¬ 
itive charge increases its maximum of 0.4 C/cm® at the 
surface at time t = 200 ps and the negative charge has 
the maximum of —0.23 C/cm® in the center of impact 
zone at time t = 18 ns. 

Obviously, the electric potential closely follows the dis¬ 
tribution of the total charge. Initially we observe a pos¬ 
itive potential in the middle of the impact zone sur¬ 
rounded by a shell of weak negative potential. For ex¬ 
ample, in Si 02 at the beginning a positive potential with 
the maximum of O.ll V is observed stretching across the 
sample-vacuum interface. A shell of weak negative po¬ 
tential is situated around this positive central region and 
appears inside the sample only. After a few hundred pi¬ 
coseconds, a transition occurs which results in a different 
situation for potential and that is the negative potential 
appears in the middle with shells of positive and negative 
potentials, respectively. The minimal value achieved by 
the potential during the first microsecond is —0.1 V in 
Si 02 which is situated in the center of impact zone and 
in the period of 20 to 50 ns. 



FIG. 5. Evolution of the total charge {p + pt — n — nt)q 
(C cm“®) in Si02 during bombardment with a beam current of 
160 nA. Top row: vertical cross-section, side length is 100 nm; 
bottom row: top-view of the sample-vacuum interface, diam¬ 
eter is 200 nm. 



FIG. 6. Evolution of the total charge {p -\- pt — n — nt)q 
(G cm“®) in Si02 during bombardment with a beam current of 
160 pA. Top row: vertical cross-section, side length is 100 nm; 
bottom row: top-view of the sample-vacuum interface, diam¬ 
eter is 200 nm. 


C. Electron bombardment 

The electron gun of a typical SEM is able to produce 
PE currents in the range of pico to nano Amperes (i.e. 
average interval between PE impacts from nano to pi¬ 
coseconds). The charge dynamics following the impact 
of a single PE, analyzed in the previous section, clearly 
shows that the next electron faces highly variable condi¬ 
tions in the sample depending on the time of its arrival. 

Since the main features of the charge dynamics in 
AI 2 O 3 and Si 02 are essentially similar, we restrict our 
discussion to the latter material. In this section a Si 02 
sample is considered under focused beams with the cur¬ 


rents of 160 nA and 160 pA (average times between PE 
impacts are 1 ps and 1 ns, respectively). To illustrate 
the nontrivial effect of the varying PE current the results 
in Figures and 1^ are presented for the same number of 
PE impacts in both beams that, obviously, correspond to 
different illumination times. 

We start with the higher current of 160 nA modelled 
as a sequence of PE’s arriving with exact one picosecond 
intervals between them. Figure shows the evolution of 
the total charge density during the first 50 impacts. On 
the fine temporal scale (not shown) we observe that the 
densities of free electrons and holes reach their maxima 
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of 9.05 X 10^® and 1.06 x 10^° cm“^, respectively, at the 
end of the generation (ionization) stage after impacts and 
decrease afterwards. 

The maximum density of trapped electrons reaches its 
maximum (the density of trapping sites iV„) for the first 
time at t = 30 ps (after 30 PE impacts) inside the impact 
zone. A similar local saturation for the trapped holes 
does not happen during the first 50 PE’s, although, their 
density grows. 

Similar to the aftermath of a single PE impact we ob¬ 
serve a (semi) spherical wave of charge density emerg¬ 
ing from the impact zone. However, now it remains a 
growing positive charge zone surrounded by the shell of 
negative charge without the charge-sign oscillation as in 
Eigure At the very beginning the positive charge has 
access to the surface, but towards the 50th impact a 
layer of negative charge prevents the positive charge from 
touching the surface. The maximum positive charge of 
1.36 G/ctc? is observed at the end and at the surface 
(lower-left image of Figure]^. The positive charge in the 
center of the expanding zone increases to 0.05 C/cm^ 
towards the end (upper-right image of Figure]^. The 
negative charge remains confined to a shell around the 
positive charge. This shell becomes distorted by grow¬ 
ing thicker with time along the sample-vacuum interface 
with the distance from the injection point, thus, reaching 
the maximum value of —0.08 G/cve? at the 50th impact. 
In the beginning and at time 8 ns, the postive potential 
reaches its maximum of 0.11 V in the center of impact 
zone and the negative potential reaches the minimum of 
—0.15 V at the end of this initial bombardment period. 

Next, we consider the beam current of 160 pA corre¬ 
sponding to one nanosecond intervals between PE im¬ 
pacts. Again, for a very short time after each impact, 
an increase in the density of free electrons is observed, 
which, after less that 0.5 ns, drops to the almost the 
intrinsic carrier density of the material. The maximum 
of 2.08 X 10^® cm“® occurs at the end of the generation 
(ionization) stage. The density of free holes reaches its 
maximum of 4.18 x 10^® cm“® in the middle of the gen¬ 
eration stage. 

For about 30 ps after each impact an increase in the 
density of the trapped electrons is seen after which the 
density remains constant until the next impact. Compar¬ 
ing this with the previously considered higher current we 
observe that now it takes a longer time (21 ns) but less 
impacts (21 PE impacts) for the density of the trapped 
electrons to reach the density of trapping sites in the mid¬ 
dle of the impact zone. Similarly to the previous higher 
current, the density of trapped holes does not reach the 
trapping site density during the considered bombardment 
period. 

Comparing the surface charge of the high (Figure 
bottom) and low (Figure bottom) currents we see that 
with the higher current, the positive charge at the in¬ 
jection point is present surrounded by a ring of negative 
charge, whereas with the lower current the charge at the 
injection point is negative surrounded by a ring of posi- 



FIG. 7. Electron flux through the sample-vacuum interface 
in Si02 illuminated by the 160 nA beam obtained with the 
pulsed model (regular and random impacts) and the time- 
uniform source model. 

tive charge (Figure [^bottom). In fact, with the 160 pA 
current the positive surface charge is also seen, but not 
at the time of impact. Subsequently, after each impact, 
the positive potential is observed at the center of impact 
zone and after a few hundred picoseconds it is replaces 
with a negative potential. In the middle of the ionization 
time, both positive and negative potentials reach thier 
maxima of 0.17 V and -0.07 V, respectively. 

Apart from the absence of rapid charge oscillation in 
the expanding spherical wave pattern, one can notice an¬ 
other substantial difference with the dynamics of Fig¬ 
ure Namely, the visible speed of expansion of the 
charged zone is much slower under the bombardment con¬ 
ditions, than during a single impact. Later we shall see 
that this speed of expansion roughly corresponds to the 
growth of the zone occupied by the trapped charges. 

We have extended the simulation to 500 PE impacts 
(at 160 nA) and compared the idealistic source with pe¬ 
riodic impacts considered above with the more realistic 
source whose PE’s impact the sample at time instants 
drawn from the Poisson distribution. Another purpose 
of this 500-impact simulation is to test the applicability 
of the time-uniform source model, see eq. (^. Figure 
shows the electron flux through the sample-vacuum in¬ 
terface obtained by the pulsed model with regular and 
random impacts as well as time-uniform source model. 
The result shows that for this relatively short period at 
the beginning of the bombardment all three fluxes are in 
good agreement with each other. Based on these promis¬ 
ing results, simulations for longer intervals of time can 
be carried out with the time-uniform source model. 


D. Steady state 

From the mathematical point view no steady-state so¬ 
lution exists with the pulsed source model (where each 
PE impact is modelled individually). At most, one can 
expect a time-periodic solution if PE impacts happen at 
regular intervals. The time-uniform source may, on the 
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FIG. 8. Electron flux through the sample-vacuum interface 
(top plot). Insert shows the zoomed-in portion of the plot cor¬ 
responding to the end of the simulated interval, where the os¬ 
cillating curve represents flux obtained with the pulsed-source 
model. Lower images: the evolution of the surface charge ob¬ 
tained by the time-uniform source model with the 160 nA 
beam current. 


n (1 ns) 




P,(1 ns) 





(400 ns) 



FIG. 9. The evolution of the densities of trapped electrons 
and holes obtained by the time-uniform source model with the 
160 nA beam current. Vertical cross-section, width - 100 nm, 
height - 200 nm 


other hand, result in a solution that converges to a steady 
state for t —oo. In this section we investigate the large 
t behaviour of the time-uniform source model in various 
circumstances. 


Figures [8 -11 show the simulation results by the time- 
uniform source model in the case of the high beam cur¬ 
rent of 160 nA. The electron flux through the sample- 
vacuum interface in Fig. shows that a steady state 
starts around 400 nanoseconds in this case. To con¬ 
firm that this, indeed, is a steady state we initiate the 
pulsed-source bombardment of the sample with the ini¬ 
tial conditions set to the time-uniform steady-state so¬ 
lution. If such a steady state is stable, then the corre¬ 
sponding charge distributions and potential could be in 
the neighbourhood of a time-periodic solution expected 
in the case of the pulsed-source model with impacts at 
regular intervals. As can be seen from the magnified por¬ 
tion of the flux plot in Fig. the flux computed with the 
pulsed-source model, indeed, appears to oscillate around 
the time-uniform flux. 


The evolution of the surface charge shown in Fig. 
starts with the positive charge at the injection point sur¬ 
rounded by a ring of negative charge. However, as the im¬ 
ages show, after a few nanoseconds, this negative charge 
ring is removed by an outward-going (surface) wave and 
the positive charge settles in the center as a steady state. 
This positive charge grows up to the value of 88 C/cm^. 

This behaviour and other features of the total charge 
dynamics are strongly influenced by the evolution of the 
trapped charge density shown in Fig. As one can 


see, the trapped charges rapidly reach their saturation 
value Nn in the impact zone, after which this satu¬ 
rated trapped-charge zone spreads outwards towards the 
ohmic boundaries. The electron trapping process and the 
spread of the corresponding zone is slightly faster than 
that of the holes (due to the higher mobility and trap¬ 
ping cross-section of the electrons). This, in particular, 
explains the presence of a slowly spreading negative shell 
in the top images of Figures [5]-{^ Also, around t = 100 ns, 
a surface channel reaching the omic contacts developes, 
consisting of saturated trapped holes. It provides a path 
free of trapping for the positive charge transport along 
the sample vacuum interface. The latter can explain the 
saturation of the positive surface charge density shown 
in Fig. [8| 

Further, the steady state with this beam current is 
characterized by the densities of free electrons and holes 
around 2.31 x 10^^ cm“^ and 2.49 x 10^^ cm“^, respec¬ 
tively. The evolution of potential is shown in Fig. |10| 
It can be noted that potential follows the surface charge 
behavior, starting as a negative ring around a positive im¬ 
pact region. Then, the negative potential moves to the 
ohmic boundaries with time, and, after a few nanosec¬ 
onds, the positive potential dominates in the sample as 
well as the vacuum close to the interface. A weak nega¬ 
tive potential is observed below the positive one, however, 
it disappears after a few nanosecond. 

Since Nn = Np in the present simulations, and both 
the trapped electrons and the trapped holes reach their 
saturation values almost everywhere (thus, effectively 
neutralizing the total trapped charge), the positive po- 
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1 ns 20 ns 400 ns 



FIG. 10. The evolution of potential obtained by the time- 
uniform source model with the 160 nA beam current. Vertical 
cross-section (vacuum part included), width - 100 nm, height 
- 400 nm 



z (/im) 

FIG. 11. Evolution of the total charge along the beam direc¬ 
tion (z;-axis) with the 160 nA beam current. The coordinate 
z; = 0.2 corresponds to the sample-vacuum interface. Vertical 
axis is clipped. See Fig. [8] for the actual surface charge values. 


tential at large t is the result of free rather than trapped 
holes. For example at t = 400 ns, the free charge, {p—n)q, 


at the center of the surface is about 89.4 G/cni while the 
trapped charge, (pt — nt)q, at the surface and around the 
injection point is almost homogeneous and around zero 
and close to the ohmic contact is 2.56 C/cm^. 

The evolution of the total charge along the beam di¬ 
rection when the beam current is 160 nA is shown in 
Fig. where the charge is displayed at three differ¬ 


ent time instants. The slight spatial advance of the sat¬ 
urated trapped-electron zone with respect to the satu¬ 
rated trapped-holes zone may explain the slowly expand¬ 
ing negative charge shell. As mentioned above, in the 
equilibrium state the sample is positively charged by the 
free holes. At this stage the trapped electrons and holes 
cancel each other, whereas, the free electrons disappear 
not only at the ohmic contacts, but through the sample- 
vacuum interface as well. 


Figure [T^ shows the simulation results with the lower 
160 pA beam current. The electron flux through the 
sample-vacuum interface demonstrates that the system 
reaches the steady state at a later time, compared to the 
previous higher 160 nA current, but with fewer PE im¬ 
pacts. A significant difference is observed between the 
free charge distributions for these two values of beam 
currents. In particular, the shapes of the free electrons 
and free holes distributions in the steady state are same 
in the case of higher 160 nA beam current. Whereas, the 
difference between these distributions at 160 pA is obvi¬ 
ous from the images of Fig. |12| The maximum density 
of free electrons is at the center of the impact zone, while 
for the free holes it is at the surface. The overall positive 
potential in the steady state is the result of the excess of 
free holes at the surface. 

The total charge distribution along the beam direction 



FIG. 12. Electron flux through the sample-vacuum inter¬ 
face, electric potential (vertical cross-section, vacuum part 
included, width - 100 nm, height - 400 nm), and the densities 
of free charges (vertical cross-section, side length - 100 nm) 
obtained by the time-uniform source model with the 160 pA 
beam current. 


for the beam current of 160 pA is shown in Fig. It is 
apparent that the spatial variations in the total charge 
density happen closer to the sample-vacuum interface (at 
depths less than 70 nm) if compared to the higher 160 nA 
beam current. Also, the surface charge remains positive 
from the start of the process. 

Finally, we consider an even lower 16 pA beam current. 
The electron flux through the sample-vacuum interface 
in Fig. and comparison with the previously obtained 
fluxes due to higher beam currents bring us to the conclu¬ 
sion that the flux rate and the time it takes to reach the 
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FIG. 13. Evolution of the total charge along the beam direc¬ 
tion (^-axis) with the 160 pA beam cnrrent. The coordinate 
z = 0.2 corresponds to the sample-vacuum interface. 


steady state are roughly proportional to the beam cur¬ 
rent. Yet the number of PE impacts to reach the steady 
state is roughly inversely proportional to the beam cur¬ 
rent. The total charge distribution along the beam di¬ 
rection shown in Fig. demonstrates that, compared 
to the 160 pA current, the signihcant spatial variations 
happen closer to the sample-vacuum interface, approxi¬ 
mately within the depth of 40 nm. The system reaches 
the steady state around 0.1 ms. The densities of free elec¬ 
trons and holes reach the maxima of 1.77 x 10^® cm“^ 
and 4.5 x 10^^ cm“^, respectively. 


This simulation confirms that a lower beam current re¬ 
sults in signihcant differences in the densities and spatial 
distribution between the free electrons and the free holes. 
At lower beam current the free electrons tend to concen¬ 
trate at the center of the impact zone at some distance 
to the sample-vacuum interface, while the free holes are 
densely concentrated at the interface (see Fig. 14). 

The lower-right image of Fig. 14 shows the spatial dis¬ 
tribution of the recombination term U, with the high¬ 
est rate of 8.25 x 10^^ cm“^s“^ achieved in the steady 
state. In the beginning of the process the recombina¬ 
tion occurs at the center of the impact zone and at 
the sample-vacuum interface, approximately at the same 
rate. After a short time (few microseconds), the recom¬ 
bination is mostly active around the center of the impact 
zone due to the low density of free electrons at the sur¬ 
face. The highest recombination rates of 5.45 x 10^® and 
5.79 X 10^® cm“®s“^ in the steady state are obtained 
with the beam currents of 160 pA and 160 nA, respec¬ 
tively. These results indicate that the recombination rate 
increases almost linearly with the beam current. Also, 
the variation of the recombination rate with the beam 
current clearly shows that applying a constant recombi¬ 
nation rate in computational models does not reflect the 
actual recombination process in dielectric samples under 
SEM. 



FIG. 14. Electron flux through the sample-vacuum interface 
and the evolution of the total charge along the beam direction 
(z-axis) with the beam current of 16 pA. The bottom row 
shows the densities of free charges and the recombination rate 
in the steady-sate regime (vertical cross-section, side length - 
100 nm). 


The SE yield can be calculated from the electron flux 
through the sample-vacuum interface shown in Figures]^ 
[T 2 I and A clear discrepancy in the measured SE 
yield of insulators is found in the literature, which can 
be attributed to the differences in the assumptions and 
conditions of the experiments. The present simulations, 
where the only varying parameter is the applied beam 
current, show that the SE yield increases with the beam 
current (provided all other conditions are fixed). The SE 
yields obtained here for the particular Si 02 sample un¬ 
der focused beam currents of 16 pA, 160 pA, and 160 nA 
are: 0.12, 1.67, and 2.74, respectively, which are, gener¬ 
ally, within the range of experimental values reported in 
the literature. In fact, we observe a weak (logarithmic) 
dependence of the SE yield on the beam current. The 
experimental values of the SE yield for the Si02 (steam 
formed) sample in the database by JojElare: 0.25, 1.02, 
and 1.18. The measured SE yield for the Si02 (quartz) 
sample iiPis approximately 3. Whereas, according to the 
experiment of Yong et al.l^, the SE yield of a “wet and 
sputtered” silicon dioxide is found to be greater than 3. 


CONCLUSIONS 


We have proposed a fully self-consistent drift-diffusion- 
reaction model augmented with a dynamic charge trap¬ 
ping model for the quantitative numerical investigation 
































14 


of the electron beam interaction with dielectric samples. 
The pulsed and non-equilibrium nature of the charge in¬ 
jection mechanism and the back reaction of the accumu¬ 
lated charge on the incoming primary electrons are incor¬ 
porated in the model through an explicit semi-empirical 
source formula. We have presented and compared two 
approaches to the charge injection problem. The first 
one is a pulsed source model reflecting the actual dis¬ 
crete nature of the electron beam. The second approach 
reduces the computational burden by applying a tem¬ 
poral average of the actual pulsed source function, which 
allows simulation at much longer time scales. Our results 
confirm the agreement between these two approaches in 
the initial stage and in the saturation regime. 

The proposed model features a Robin-type semi- 
permeable boundary condition at the sample-vacuum in¬ 
terface reflecting the fact that the electrons are allowed 
to go through the boundary, while holes are not. We 
have devised a simple optimization procedure to deduce 
the surface recombination velocity (SRV) of dielectrics in 
vacuum from the experimental SE yield data. 

The results of our simulations clearly demonstrate the 
need for the dynamic trapping equations in computa¬ 
tional models of this kind. The trapping dynamics, 
namely, the time evolution of the spatial distributions 
of the trapped charge densities, has a major influence 
on (and helps to explain) the total charge distribution 
within the sample and the apparent transients in the sur¬ 
face charge density (in the high-current regime). 

Inclusion of the dynamic generation-recombination 
model is also deemed necessary, since, as it turns out, 
the recombination rate depends on the beam current. 
Other quantities depend either strongly or weakly on the 
beam current as well, e.g., the local charge densities in 
the steady state show a linear dependence, whereas, the 
dependence of the SE yield turns out to be logarithmic. 

Another conclusion of our study that requires a deeper 
mathematical analysis is the apparent existence of a time- 
periodic steady state in the considered system of equa¬ 
tions for a pulsed source with PE impacts at regular in¬ 
tervals observed in the neighborhood of the steady-state 
solution for an averaged time-uniform source. 


ACKNOWLEDGEMENTS 


The authors acknowledge the partial financial support 
of the FEI Company (Netherlands) and many fruitful dis¬ 
cussions with Dr. S. Sluyterman and Dr. E. Bosch (both 
with FEI). Thanks are also due to Professor B.J. Thijsse 
(Delft University of Technology), who gave us valuable 
counsel in the early stages of this work. Finally, we thank 
the anonymous reviewers for many thoughtful comments 
and suggestions. 


Appendix A: Weak form 


Consider the partial differential equations JB, 0, and 
@ in the axisymmetric geometry (see Fig: 
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To derive the weak formulation we integrate over the 
cross-sectional area (rdrdz) arriving at the following form 
of the equation (Al): 
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where w is the weight function. Integrating the highest- 
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where v = is the outward unit normal vector to 

the boundary. Therefore, the weak form of the equation 
(Al) can be written as: 
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or 


Substitution of the boundary conditions gives: 
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The weak form of the equation (A2) is: 
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Integrating the highest-order terms by parts we get: 
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Therefore, the weak form of the equation (A2) can be 
written as: 
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Along similar lines the weak form of the equation (A3) 
cab be derived as: 
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where p |an 2 = and w |an 2 = 0. 
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